#include <iostream>
#include <vector>
#include <cmath>
#include <sstream>
#include "../Eigen/Dense"
#include "../spline_functions.h"
// #include "matplotlibcpp.h"

// namespace plt = matplotlibcpp;

using namespace std;
using namespace Eigen;

// 定义函数 f(x)
double f(double x) {
    return 1.0 / (1 + 25 * x * x);
}

// 计算插值误差的最大范数
double max_norm_error(const vector<double>& x, const vector<double>& y, int N, double m1, double mN) {
    vector<double> y_spline(N);
    ClampedCubicSpline_p spline = ClampedCubicSpline_p(x, y, m1, mN);

    double max_error = 0.0;
    double h = 2.0 / (N - 1); // 子区间宽度
    for (int i = 0; i < N - 1; ++i) {
        double mid_point = -1 + (i + 0.5) * h; // 子区间中点
        double f_exact = f(mid_point);
        double f_spline = spline.evaluate(mid_point);
        double error = fabs(f_exact - f_spline);
        if (error > max_error) {
            max_error = error;
        }
    }
    return max_error;
}

// 绘制函数和样条插值
void plot_spline(const vector<double>& x, const vector<double>& y, int N, double m1, double mN) {
    ClampedCubicSpline_p spline = ClampedCubicSpline_p(x, y, m1, mN);
    vector<double> x_values, y_values, y_exacts;

    for (double xi = -1.0; xi <= 1.0; xi += 0.01) {
        x_values.push_back(xi);
        y_values.push_back(spline.evaluate(xi));
        y_exacts.push_back(f(xi));
    }
    // plt::named_plot("Spline Interpolation",x_values, y_values, "r-");
    // plt::named_plot("Exact Function", x_values, y_exacts, "b--");
    // plt::xlabel("x");
    // plt::ylabel("f(x)");
    // plt::legend();
    // plt::grid(true);
}

int main() {
    std::cout << "-------Problem A: Applying cubic-spline interpolation-------" << std::endl;
    vector<int> N_values = {6, 11, 21, 41, 81};
    vector<double> x, y;

    for (int N : N_values) {
        x.clear();
        y.clear();
        double h = 2.0 / (N - 1);

        for (int i = 0; i < N-1; ++i) {
            double xi = -1 + i * h;
            x.push_back(xi);
            y.push_back(f(xi));
        }
        x.push_back(1);
        y.push_back(f(1));
        double m1 = 50.0/(26*26);
        double mN = -m1;
        // stringstream title_stream;
        // title_stream << "Cubic Spline Interpolation vs Exact Function (N = " << N << ")";
        // plt::title(title_stream.str());
        // plot_spline(x, y, N, m1, mN);
        // plt::save("F:/Codefield/numerical_analysis/Programming/Chapter3/figure/PA_" + std::to_string(N) + ".png");
        // plt::clf();

        double error = max_norm_error(x, y, N, m1, mN);
        cout << "N = " << N << ", Max-norm error = " << error << endl;
    }

    return 0;
}